Possible discovery of a nonlinear tail and second-order quasinormal modes 

in black hole ringdown 
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We investigate the nonlinear evolution of black hole ringdown in the framework of higher-order 
metric perturbation theory. By solving the initial-value problem of a simplified nonlinear field model 
analytically as well as numerically, we find that (i) second-order quasinormal modes (QNMs) are 
indeed excited at frequencies different from those of first-order QNMs, as predicted recently. We 
also find serendipitously that (ii) the evolution is dominated by a new type of power-law tail at 
late times. This "second-order power-law tail" decays more slowly than any late-time tails known 
in the first-order (i.e., linear) perturbation theory, and is generated at the wavefront of the first- 
order perturbation by an essentially nonlinear mechanism. These nonlinear components should be 
particularly significant for binary black hole coalescences, and could open a new precision science in 
gravitational wave studies. 

PACS numbers: 04.30.Tv, 04.25.Nx, 04.70.-s 



I. INTRODUCTION 

Direct detection of gravitational waves is one of the 
most exciting challenges in astrophysics today. It will 
not only enable to verify general relativity, but will also 
provide a new observational window toward the universe. 
Current and future gravitational wave detectors such as 
LIGO, LISA, and DECIGO/BBO should make it a real- 
ity in the near future. 

The "ringdown" of black holes is an important tar- 
get of gravitational- wave observations. It is known the- 
oretically that black holes have characteristic oscillation 
modes, called quasinormal modes (QNMs) |2|. QNMs 
are the solutions to perturbation equations which sat- 
isfy so-called "outgoing boundary conditions," i.e., which 
propagate purely inward at the horizon and purely out- 
ward at spatial infinity. They are distinguished from 
ordinary normal modes because they decay at certain 
rates, having complex frequencies. The remarkable prop- 
erty of the black hole QNMs is that their frequencies are 
uniquely determined by the mass and spin of black holes. 
This means that the parameters of black holes can be di- 
rectly measured by observing their ringdown through the 
gravitational waves. 

The black hole ringdown is particularly prominent in 
binary black hole mergers Recently, numerical rela- 
tivity has succeeded in calculating their entire evolution 
[H HI @|- As a result, it has been found that the least- 
damped QNM dominates the gravitational waves radi- 
ated during the ringdown phase, carrying away ~ 1% 
of the initial rest mass energy of the binary system [?J . 
The ringdown wave with such large energy will be de- 
tected by future gravitational-wave detectors with a high 
signal-to-noise ratio 3 . 

Usually, the black hole ringdown is understood as 
the response of a black hole to linear perturbations 
d, [T(| HH, EH- In black hole perturbation theory, 
the metric is divided into a stationary black hole back- 



ground (formed after a binary merger, for instance) and 
small perturbations on it. Now let us consider an initial- 
value problem 0] where the perturbations are initially 
distributed near the horizon. Then, the linear theory 
predicts that an observer far away from the black hole 
sees two different stages in the evolution of gravitational 
waves. At early times, the evolution is dominated by 
quasinormal (QN) ringing, a superposition of QNMs. At 
late times, in contrast, the evolution is characterized by 
a power-law tail, which is related to the backscattering of 
the QN ringing far away from the black hole [3. Also, 
the initial-value analysis clarifies another important fact 
that the QNMs are excited only around a black hole [l^] • 
As a result, evolved QNMs are always truncated by a 
future-directed light cone [l4|. This truncation is essen- 
tial for QNMs, since without it they would diverge at 
both spatial infinity and the horizon. 

Although the linear theory is clearly useful, the re- 
cent numerical simulations imply that the nonlinearity 
of the ringdown is also worth studying. The radiated 
energy of ~ 1% x (rest mass energy) is translated into 
the amplitude ~ 10% of the first-order (linear) metric 
perturbations, implying that the next-leading, second- 
order ones will have amplitude as large as ~ 10% of the 
first-order amplitude. This strongly suggests that the 
second- and higher-order perturbations should contribute 
the ringdown waveform to that extent, providing us with 
some additional information on the merger events and 
the nonlinearity of general relativity. 

Despite its potential importance, the nonlinearity of 
black hole ringdown has been hardly investigated. Only 
recently, one of the authors (K. I.) and his coworker have 
for the first time applied second-order metric perturba- 
tion theory to studying QNMs of Schwarzschild black 
holes [TU, [l6[. Because of the nonlinearity of the Ein- 
stein equation, two normal modes in first-order pertur- 
bations [say cos(f2ii) and cos(f22*)] couple each other, 
yielding their "sum tone" oc cos[(f2i + f^)^] and "differ- 
ence tone" oc cos[(f2i — f^H] in second-order ones. The 
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authors found that the sum tones and difference tones 
of the first-order QNMs also satisfy the outgoing bound- 
ary conditions, meaning that they are the QNMs in the 
second-order perturbation theory. They concluded that if 
these "second-order QNMs" are excited in a binary black 
hole merger, they will have enough energy to be detected 
by future gravitational-wave detectors. Detection of the 
second-order QNMs has many possible applications, such 
as measuring the distance from the source system, test- 
ing the nonlinearity of general relativity, and rejecting 
fake events in ringdown searches. 

Questions have remained unanswered. First, are the 
second-order QNMs really excited even if the first-order 
ones are spatially truncated? The previous work [TEl . [l6j 
was done in the frequency domain, where the first-order 
QNMs (or the sources of the second-order perturbations) 
are implicitly assumed to be extended from the hori- 
zon to spatial infinity. This assumption is actually in- 
valid, because the sources would diverge far away from 
the black hole. The first-order QNMs are actually trun- 
cated at a light cone, and thus the sources of the second- 
order QNMs are distributed in a time- dependent way. 
Second, does the second-order theory provide only the 
second-order QNMs? In general, first-order QNMs are 
accompanied by a power-law tail, as naturally shown us- 
ing the time-domain Green's function of the perturba- 
tions [1, [H, [HI]. On the other hand, it has been open 
whether the second-order QNMs are accompanied by any 
other wave components, especially of nonlinear origin. 
To address these questions, one must move away from 
the frequency-domain analysis and treat the evolution of 
the second-order perturbations as a time- dependent prob- 
lem. 

In this study, we give answers to the above questions 
by solving an initial-value problem for nonlinear black 
hole ringdown. Our analytic calculation is based on a 
time-domain analysis developed for linear black hole per- 
turbations [1, [l2j . As far as we know, this analysis has 
been never applied to nonlinear black hole perturbations. 
Also, we employ a simple nonlinear field model having 
key properties in common with black hole perturbations. 
Our nonlinear model is attractive, since the original equa- 
tions for black hole perturbations are too complicated to 
be treated analytically, and since direct simulations us- 
ing numerical relativity would not be accurate enough 
to resolve nonlinear evolution in ringdown phases at the 
present time [l5| . Our model enables us not only to solve 
the problem analytically, but also to confirm the results 
with fully nonlinear simulations in a good accuracy. 

Our answers to the above questions are as follows. 
First, the second-order QNMs are indeed excited by trun- 
cated first-order QNMs. Second, the second-order QNMs 
are accompanied by a power-law tail, but it is essen- 
tially different from the tails in the linear perturbation 
theory [13} . This tail, to which we shall refer as "second- 
order power-law tail," decays more slowly than any of the 
first-order tails, and even dominates the fully nonlinear 
evolution of the model field at late times. This means a 



surprising fact that the first-order theory fails to predict 
the late-time behavior of the ringdown. We also find that 
the behavior of the second-order tail is completely deter- 
mined by the spatial dependence of an effective "source" 
term for second-order perturbations. Since the effective 
source term of our model behaves in essentially the same 
manner as that of black hole perturbations, we strongly 
expect that the second-order power-law tail should ap- 
pear in real black hole ringdown, and should play a lead- 
ing role in the late-time evolution. 

This paper is organized as follows. In Sec. II, we in- 
troduce a model system for black hole perturbations and 
derive the perturbation equations used in the subsequent 
sections. In Sec. Ill, we analytically solve an initial- 
value problem of the perturbations with an approximated 
Green's function, and decompose the evolved second- 
order perturbation into different components. In Sec. IV, 
we confirm the validity of the analytic calculations by 
numerical simulations. In Sec. V, we present our conclu- 
sion and discuss several possible applications to future 
gravitational-wave studies. We also discuss briefly a re- 
lation between our second-order power-law tail and non- 
linear tails recently found in other systems. In Appendix 
A, we summarize how to construct the time-dependent 
Green's function using Andersson's "asymptotic approx- 
imation" [12]. 

II. NONLINEAR FIELD MODEL 

Before introducing our nonlinear field model, we briefly 
review how the equations of black hole perturbations are 
derived from the Einstein equation. The Einstein equa- 
tion in vacuum reads 

G[g] = 0, (1) 

where G is the Einstein tensor and g is a metric tensor. 
In metric perturbation theory, the metric g is expanded 
into a stationary background part g^ and perturbation 
parts (j = 1,2,---): 

g = g^ +eh^ +e 2 h^ + ■■■ , (2) 

where e is an expansion parameter. The wave equations 
for h^) can be obtained by expanding the Einstein equa- 
tion (TT]) perturbatively with e: 

G (1) [h (1) ] = 0, (3) 
G (i)[ h P) ]= _ G {2) [h W ift « ]i (4) 



where G^ is the linearized Einstein tensor, and 
G^ [a, b] is a nonlinear part of G and is quadratic in 
a and b. On the Schwarzschild background, the per- 
turbations h^' can be expanded into tensor harmonic 
components n£i, and Eqs. ([3]) and ^ can be reduced to 



3 



(l+l)-dimensional master equations [171. Ha. UM, 120. W\ 



(5) 
(6) 



where i/j 



(i) 



and 



tions for h 



Xtm are , respectively, the master func- 
^ and (for their definitions, see dEHH); 
r* is the tortoise coordinate [ItJ , covering from the hori- 
zon (r„ — > — oo) to spatial infinity (r* — » oo); 14 (r) is 
an effective potential term that has a peak outside the 
horizon and vanishes as r* — > ±oo (Regge- Wheeler [ItJ 

is an effective source 
,h {1] ] in Eq. flU) and 



or Zerilli [la] potential); and 5^ m 
term, originating from 



,-(2) 



essentially quadratic in 1/7™ ■ Note that the left-hand 
sides of Eqs. (0 and ^ are described by the same linear 
differential operator — d 2 + d 2 t — Vf, because they both 
originate from the linearized Einstein tensor G^ . 

What we want to do in this study is to clarify what 
types of wave components are created as the second-order 
perturbations during black hole ringdown. To decompose 
the perturbations into different types of components, we 
have to solve the evolution as analytically as possible. 
However, it would be challenging to do this by a di- 
rect use of Eq. 



many linearly independent components, each of which 
depends on r* in a complicated way (see, e.g., Eq. (13) 
in Ref. 15|). In addition, it would be also challenging to 
confirm the results with full-order computations of the 
Einstein equation ([!}, since current numerical relativity 
may not be accurate enough to resolve the evolution of 
the second-order perturbations [H, EH- Therefore, we 
prefer for our purpose to consider a simple model instead 
of the real perturbations. 

Now we introduce our model. To avoid oversimplifi- 
cation, we require the model to satisfy that (i) its first- 
and second-order perturbations obey (l+l)-dimensional 
wave equations with a potential barrier, like Eqs. ^ and 
([6]) , and that (ii) the effective source term for the second- 
order perturbation is quadratic in the first-order one. As 
one of the simplest models, we adopt a nonlinear, (1+1)- 
dimensional scalar field obeying 



d 2 +d 2 x - V{x) *(jb,x) = F(xMt,x) 2 , (7) 



where V(x) is a potential term and F(x) is a function 
specified later. We assume that V(x) has a peak at x = 
and vanishes at |x| 3> 1, as schematically shown in Fig. [I] 
Since x corresponds to the tortoise coordinate r* for black 
holes, we refer to the limits x — > —00 and x — > +00 as 
"horizon" and "infinity," respectively. We also assume 
an observer's position near infinity, x 3> 1. 

It is easily checked that this model meets the above re- 
quirements. We expand $(t, x) in terms of an expansion 
parameter e to write 
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FIG. 1: Schematic illustration of the potential barrier V(x) 
and the initial data <!>(£ = 0,x). The potential barrier van- 
ishes at \x\ 3> 1, and the initial data has a compact support 
near the potential peak located at x = 0. 



Substituting this into Eq. (JT]), we obtain an infinite set 
of perturbation equations, the first two of which read 



-d 2 t +d 2 x - V(x) 
-d\ + d 2 x - V(x) 



W =0, 

^ =S^(t,x), 



since the source term ST; contains where 



SW(t,x) =F(x) 1 



(9) 
(10) 

(11) 



is an effective source term for the second-order perturba- 
tion cj)( 2 \ It is clear from Eqs. (|5))- lfTTj) that this model 
does satisfy the requirements (i) and (ii). 

In this study, we consider an initial-value problem of 
Eq. ([7]) with the following initial conditions: 



$(< = 0,3!) =/(»), 

d t $(t = 0,x)=g(x) 



(12) 
(13) 



where f{x) and g(x) are arbitrary functions vanishing far 
away from the potential peak (\x\ 3>1), as schematically 
illustrated in Fig.[T] To solve the problem perturbatively, 
we have to divide these initial data among the perturba- 
tion quantities . Since there is no unique choice of 
dividing the data, we adopt the simplest choice such that 
the first-order perturbation theory is exact at the initial 
time: 



<pV(t = 0,x)=f(x), 
d t <f>W(t = 0,x)=g{x), 
<j) {k) {t = 0,x) = 0, 

d t (/> w (t = o,x) = 0, 



(14) 
(15) 
(16) 
(17) 



$(M) = £0 (1) (M) +e 2 (f> i2) (t,x) 



(8) 



where k > 2. Note that this choice does not guarantee 
that the first-order theory remains to be a good approx- 
imation to the full-order theory at all times. In fact, 
as we will see later, </>( 2 ) does surpass in amplitude, 
dominating the late-time behavior of $. 
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III. ANALYTIC CALCULATIONS 



and the QNM part Gq are obtained as (see Appendix A) 



Since all the perturbations <f>W (J = 1, 2, • • • ) obey lin- 
ear equations with the same operator — df + d 2 — V(x), 
their evolution is determined by a single Green's func- 
tion. The retarded Green's function G(t;x, x') for (fP> 
is the solution to 

d?+dl-V(x)\G{t-t';x,x') = 5{t-t')6(x-x'), (18) 

under the retarded boundary condition 

G(t; x, x') = 0, t<0. (19) 

Using this Green's function, the solutions to Eqs. (|9|) and 
(fTU)l with the initial conditions (jT^J) — (jTTJ) are written as 



</> (1) (t, x)= dt'dx'G(t - t'\ x, x)I{t', x'), (20) 

J J — OO 

<j>W{t,x) = 1 7°° dt'dx'G(t-t';x,x')S^(t',x'), (21) 



where I{t,x) = -f(x)S'(t)-g(x)S(t). 

If the background were "flat" (V = 0), Eq. (TTSJ) would 
be easily solved to give 



G(t-t';x,x') 



-G(t-t'-\c 



(22) 



where Q(z) is the step function. Equation (|2"2"|) means 
that an observer at (t, x) receives waves generated inside 
a past-directed light cone, t' < t — \x — x'\. Also, the 
"flatness" of the Green's function for t' < t— \x — x'\ rep- 
resents that the waves would propagate without being 
reflected or deformed. In our case, however, the back- 
ground is "curved" (V ^ 0), and some of the waves are 
reflected or deformed by the potential barrier. For this 
reason, the Green's function is much more complicated 
than Eq. (HSJ) and is conventionally divided into three 
parts [1, GJ]: 

G(t; x, x') — Gf(t; x, x') + Gq(t; x, x') + Gb(t; x, x'). 

(23) 

The "flat part" Gf represents the direct propagation of 
waves in asymptotically flat regions (\x'\ ^ 1). The 
"QNM part" Gq is related to the excitation of first- 
order QNMs. The "branch-cut part" Qb describes the 
backscattering of waves by the long-range part of the po- 
tential barrier. 

Unfortunately, the exact expression of the Green's 
function (|23p is unknown for most cases including black 
hole perturbations. Hence, we construct an approximate 
Green's function with Andersson's "asymptotic approxi- 
mation" [12| | . Using this approximation and assuming an 
observer's position at infinity (x 3> 1), the flat part Gf 



G f {t;x,x') w -~9(t- \x - x'\)Q(x + x' - t), (24) 



Gq{t;x,x') 



—Q(t-x 
2 v 



E 



X 
(t-x) 



y n (x') + c.c. 



(25) 



where {y n }^a are the first-order QNMs defined as 
the homogeneous solutions to the Laplace transform of 
Eq. (Tig)) with the "outgoing boundary conditions" 



y n {x) 



B L (s n )e~ 



x 

X 



-oo, 
-co. 



(26) 



Note that "QN frequencies" cu n , often used in the liter- 
ature, are related to s n by uj n = is n . It is also noted 
that the branch-cut part Gb is not derived from the 
asymptotic approximation, because this approximation 
implicitly neglects long-range components of the poten- 
tial barrier [121] . However, the contributions of Gb (i-e., 
late-time tails in the first-order theory) is usually small 
compared to the other contributions (see, e.g., [13!]), and 
hence we neglect Gb in the following analysis. 

The supports oiGF{t—t' '; x, x') and GQ(t — t'; x, x') are 
indicated in Fig. [U(a). For t' > t — x — \x'\, the Green's 
function G ~ Gf behaves in essentially the same way 
as that of the flat background, or Eq. (|22p. This means 
that waves generated in this region arrive at (i, x) with- 
out passing through the potential barrier, and thus with- 
out being reflected or deformed. For t' < t — x — \x'\, on 
the other hand, the Green's function G « Gq behaves 
differently from Eq. ([2"2"]) , oscillating with QN frequen- 
cies. This implies that waves generated in this region 
pass through the potential, and are then deformed into 
the first-order QNMs. These facts also imply that the 
first-order QNMs are excited only around the potential 
peak, and are thus always truncated by a future-directed 
light cone [Ti |. 

Now we perform the integration in Eq. (|20[) using 
the approximate Green's function, or Eqs. ([23]) - ([25]) and 
Gb ~ 0. Given that f(x) and g(x) are localized near the 
potential peak, the evolution of is calculated to be 



W{t,x)n-Q{t-x) J dx'^fix^dt+gix^GQ^x^') 

OO 

= Q(t - x) J2 [C n e s " {t ~ x) + C* n e< (t ~ x) 



n=0 



where 



G 



dx' 



f(x')s n +g(x') 



yn(x') 



(27) 



(28) 



is the "excitation coefficient" for the first-order QNM 



Sn(t — x) 



Equation (|27p means that the initial first-order 



5 




x' 



FIG. 2: (a) Space-time diagram describing the support of the Green's function G(t — t';x,x') [Eq. (I23|l ] for an observer 
at x (< t). supp(Gf) and supp(Gq) represent the supports of the fiat part Gf [Eq. f|24J>] and QNM part Gq [Eq. (|25l) ]. 
respectively, u' = t' — x' and v' = if + are the retarded and the advanced null coordinates, (b) Support of the source 
term 5 (2) (i',a;') [Eq. fl30}], which is equal to that of the evolved first-order perturbation w (t',x') [Eq. l[2"7))]. The stripes 
indicate the wave pattern of schematically, (c) Integration regions for 4$ and ^ in Eq. J29J, supp(G F ) n supp(5' (2) ) 
and supp(Gq) (~1 supp(5 ,( - 2 '), respectively. 



perturbation evolves into a superposition of the first- 
order QNMs truncated at x = t. Performing a similar 
calculation for x < shows that the first-order QNMs 
are also truncated at x — —t. The evolution of <f)^ is 
schematically illustrated in Fig^b). 

The evolution of (p^ is obtained by substituting 
Eq. (H71) into Eq. (JTTJ) . Since </> (1) vanishes at x > t 
and x < — £, we can rewrite Eq. (f2Tj) as 



</> (2) M)«- 



du' 



du' 



dv' G F (t-t';x,x')S^(t',x') 



dv' G Q (t-t';x,x')S^(t',x'), 

(29) 



where we have introduced the retarded and advanced null 
coordinates, u = t — x and v = t + x, and used the 
relation dt'dx' = (l/2)du'dv'. <j${t,x) and 4> { Q{t,x) 
denote the contributions from the flat part and QNM 
part, respectively. Their integration regions at t > x are 
illustrated in Fig. [2jc) . 

To carry out the integrations in Eq. ([2"9")) , we must spec- 
ify the form of the source term S^- 2 \t, x), or F(x). Since 
the (regularized) source term for black hole perturba- 
tions [Sf*](t. r») iu Eq. (G)] behaves as oc r~ 2 at infinity 
[HI flr]| . we choose F(x) so that S^(t, x) behaves as 
oc x~ 2 at x — -> oo. As one of the simplest form, we as- 
sume F(x) = (\x\ + 1)~ 2 , or 



SW(t,x) = 



^{t.xf 



(N + i) 2 ' 

which behaves as oc x~ 2 at x — > ±oo f22| . Substituting 
Eqs. ^2TJ) and (J30J) into Eq. ((29j) , we have 4> {2) {t,x) = 
for u < 0, and 

<j}f(t,x) « -- / dw' / du' 



n.n'— 



(x' + l) 2 



(31) 



for u > 0. This integration is formally performed by 
introducing an analytic function 



T{z) = e z Ei(-z), 



(32) 



where 



Ei(z) = — I dw 

w 



(33) 



is the exponential integral. Using T(z), we find a useful 
equality 



U pV 



1 r du' r dv' ' 



4./o " J u " [x' + lf 



\t{s{v -u' + 2)) - T(s(u -u' + 2)) } 



«'=0 

(34) 



(30) 

Hence, we can rewrite Eq. (j3"Tj) as 



<j>f (t, x) « 0^ (f , z) + 4 2J (t , z) , (35) 



(2) (2) 

where 4> 2 q and 0^ are the contributions from the inte- 
rior (u' = u) and the edge (u' = 0) of the source distri- 
bution [see Fig. [21(c)] : 
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<${t,x) =~J2 C n C n ,e( s " +s n'W '{t((s„ + s n ,)(v - v! + 2)) - T((s n + s n ,)(u - v! + 2))} + 

n,n' 

= - C„C„'e (s " +s "' )u {T 1 (2( S „ + s n ,)(x + 1)) - T(2(s n + s„0)} + • • 



(36) 



4 2) (*. *) = + C n C n ,e^+ s ^ u ' [T((s n + s n ,)(v -u' + 2)) - T((s„ + a n ,)(« - u' + 2))} + 

n.n' 

)(« + 2))-T(( )(« + 2))} + -. 



u'=0 



(37) 



where "■ • ■ " represents the terms proportional to C n C*, 
or C*C*/ . The behavior of these terms is essentially the 
same as that of the C n C n * terms. 

(2) (2) 

We are particularly interested in how 4>\q and (fr T 
behave at the observer's position, or at a; > 1. First, 

(2) 

we examine the behavior of 4> 2 q- Using the asymptotic 
expansion of T(z) at \z\ ^> 1, 



(38) 



fc=i 



we have 



C nn ,e^ +S ^ u + D nn ,e^+<'^ + c.c. 



(39) 

for x 3> 1. Here we have defined C nn i and D nn i as 

C nn ,=C n C n ,T(2(s n + s n ,)), (40) 

£) nn ^C n C:,T(2( Sn + <,)). (41) 

and neglected the terms of 0(x~ r ). It is helpful to rewrite 
Eq. (3gj as 



43( 



+ D„„,e- 4(0 "-°"'^ + c.c.} 



(42) 



where the 7„ = Re(s n ) < and £l n = |Im(s n )| are the 
damping rates and central frequencies of the first-order 
QNMs, respectively. Clearly, the first and second terms 
in square brackets are the "sum tone" and "difference 

(2) 

tone" of the first-order QNMs, meaning that 4> 2 q is a su- 
perposition of the second-order QNMs. That is, Eq. (j4"2"|) 
proves that the second-order QNMs are indeed excited 
during the evolution of the nonlinear field. It is also clear 
that C nn i and D nn i represent the excitation coefficients 
of the second-order QNMs. 



(2') 

Next, we examine the behavior of <j) T . Using Eq. I|38[). 
we find that this term approximately behaves as a power- 
law tail, 



^\t,x) « H 



1 



1 



t + x + 2 



where 



Sri i &n f 



+ c.c. 



(43) 



(44) 



is the excitation amplitude of the tail. Equation 
implies that this tail behaves as i -1 for x < t <, 2x, and 
as t~ 2 for t <; 2x (note that we have assumed that x ^> 1, 
but not assumed that f>x). 

Interestingly, this power-law tail is essentially different 
from the tails known in first-order perturbation theory 
(l3| . First, the tails in the first-order theory originate 
from the branch-cut part Gb of the Green's function 
(which has been neglected in our analysis), whereas our 
power-law tail appears from the flat part Gf- Second, the 
first-order tails decay more quickly than t~ 2 [Hj], while 
our tail decays more slowly than t~ 2 . Third, the behavior 
of the first-order tails depends on the form of V(x), but 
that of our tail is independent of it. In fact, as we will 
see later, it is determined by the long-range behavior of 
the source term S^(x), or F(x). This fact implies that 
our tail is generated by an essentially nonlinear mecha- 
nism. To distinguish our tail from those in the first-order 
theory, we call it the second-order power-law tail. 

In the above analysis, we have not computed the con- 
tribution from the QNM part of the Green's function, or 

(2) 

4>q ■ One would want to evaluate this by substituting the 
approximate form of Gq [Eq. ([25]) ] into Eq. (f2"Tj) . Unfor- 
tunately, such a calculation would not lead to a correct 
answer. The problem is that the dominant contribution 

(2) 

to 4>q comes from x ~ 0, where the asymptotic approxi- 
mation becomes unreliable. Nevertheless, we can expect, 

(2) 

from the physical meaning of Gq , that (fig will be a lin- 
ear combination of the first-order QNMs (oc e s "(*^ ;r )). In 
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fact, as shown in the next section, the first-order QNMs 
appear not only in (f>^ , but also in </>( 2 ) . It is noted, how- 
ever, that first-order QNMs appearing in 0( 2 ' are not im- 
portant, since they behave the same as the first-order per- 
turbation (f)^. On the contrary, the second-order QNMs 

(2) . (2) 

(f> 2 Q and power-law tail (jr T are important, because they 
behave differently from <j>^ . 

It will be helpful to mention the origins of the second- 
order QNMs and tail. In the first-order perturbation 
theory, QNMs are produced by the QNM part Gq of 
the Green's function, and power-law tails by the branch- 
cut part Gb- This means that they are both gener- 
ated in curved regions: the first-order QNMs are excited 
around the peak of the potential barrier and the first- 
order tails are yielded by the long-range decay of the po- 
tential. On the other hand, the second-order QNMs and 
power-law tail arise from the "flat part" Gf- This means 
that they are generated in asymptotically flat regions. In 
particular, the second-order power-law tail comes from 
the "edge" of the source distribution, or the truncation 
of the first-order QNMs, at u = [see Eq. and 
Fig. H£c)]. Therefore, the second-order tail cannot ap- 
pear in frequency-domain analysis, where the truncation 
is not taken into account. 

With this fact, we can easily understand why a 
nonoscillating tail can be produced by oscillating QNMs. 
As illustrated in Fig. EJc), an observer at (t,x) "sees" 
the edge extending from (V, v') = (0, t — x) to (u' , v') — 
(0, t + x) (the observer could not "clearly" see the edge 
at v' < t — x, because waves generated in this region 
are partly scattered by the potential barrier before arriv- 
ing at the observer). Since the first-order perturbation 
is constant along this edge, waves generated there also 
have a constant phase and hence do not oscillate. 

It is also easy to see that the form of the second-order 
tail is completely determined by the spatial dependence 
of the source term, i.e., F(x). In fact, the second-order 
perturbation generated at the edge v! = can be simply 
calculated to be 



~J du ' J t + dv'(t>W(t',x') 2 F(x')6(u') 

nt+X 1 1 

oc / dv'Ftv' 72) oc -, 

J t _ x x ' ' t-x + 2 t + x + 2 



(45) 



where we have used that = constant and x' = t' = 
v'/2 along the edge, and that F(x') = (\x'\ + 1)~ 2 . It is 
clear from Eq. (|4"5")1 that the second-order power-law tail 
in Eq. (|43|) is essentially the integral of F. 

We summarize the main results of our analytic calcu- 
lations. We have proven that the second-order QNMs do 
appear in the solution to the second-order perturbation. 
Also, we have found that the excitation of the second- 
order QNMs is accompanied by a slowly decaying power- 
law tail. This "second-order power-law tail" is essentially 
different from those known in the first-order theory. In 
particular, its behavior is determined by the source term, 



not by the potential term. This means that the second- 
order tail is of nonlinear origin. 



IV. NUMERICAL SIMULATIONS 

The above analysis has been largely dependent on two 
approximations: the perturbative approximation up to 
second-order and the asymptotic approximation used to 
construct the Green's function. In this section, we nu- 
merically confirm that these approximations are actually 
valid for an appropriate initial condition. 



A. Setup and method 

We performed two types of simulations. In the first 
type, we integrated the first- and second-order pertur- 
bation equations, (HJ) and (fT0|) . In the second type, we 
integrate the full-order equation ([7]). For both types of 
simulations, we use the same initial condition, potential 
term, and source term. 

As an initial condition, we chose a momentarily sta- 
tionary Gaussian wave packet 



f{x) = exp[-(2.0x) 2 



g(x) = 0. 



(46) 



As the potential term V(x), we adopted the Poschl- 
Teller potential 



Vpt(x) 



V 



cosh 2 (Kx) ' 



(47) 



where Vo and K~ x are the height and the width of the po- 
tential barrier, respectively. It is noted that the Poschl- 
Teller potential damps exponentially at both x — ► ±oo, 
while the effective potential for Schwarzschild black holes 
[Vi(r*) in Eqs. © and ©] decays as r~ 2 at — » +oo. 
However, this difference does not affect the essential fea- 
tures of the second-order QNMs and power-law tail, since 
the flat part Gf of the Green's function remains un- 
changed as long as the potential asymptotically vanishes 
far away from its peak ([H|, [23[; see also Appendix A). 

The QN frequencies for Vpt are analytically known to 
be |3 

/ 1\ / K 2 \ 1 / 2 
s n = -K[n+-)-i{y Q -—) . (48) 

In the simulations, we chose the values (Vq,K) = 
(5.0, 1.0); the "frequency" s of the least-damped (n = 0) 
first-order QNM is calculated from Eq. ([35]) to be 



so 



-0.50-2.18i, 



(49) 



or, in a more familiar notation, ujq = isq = 2.18 — 0.50z. It 
follows that the frequencies of the least-damped second- 
order QNMs are so + sq = 2sq = 2-fo — 2iflo for the "sum 
tone" mode, and so + Sq = 270 for the "difference tone" 
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mode. We use these values to compare the numerical 
results with the analytic results obtained in the last sec- 
tion. We do not use the overtone (n > 1) modes, since 
their contributions are usually smaller than that of the 
least-damped modes and since comparison would become 
much more difficult. 

In integrating the wave equations, we used a stan- 
dard second-order finite difference scheme. The com- 
putational domain was taken to be < t < 50.0 and 
—60.0 < x < 60.0, so that the wave fronts of the evolved 
perturbations never reach the spatial boundaries. The 
mesh size (At, Ax) was set to (0.01,0.1), so that the 
Courant number v = At/ Ax satisfies the stability condi- 
tion v < 1. 



B. Result 1: The first-order and second-order 
perturbations 



In Fig. 1 we show the numerical waveform $^fo s (t) = 
(t, a; bs) of the first-order perturbation observed at 
x = x hs = 5.0. We note that its wave front arrives at 
t « 4.0 ~ x bs — 1-0, not just at t — x b s . This is because 

the initial wave had a width of Sx ~ 1. From the result 

(i) 




of Sec. Ill, </> obs (t) should be a superposition of the first- QNMs. To decompose ^ s (t) into these components, 



FIG. 3: Comparison of the numerical waveform of the 
first-order perturbation, 4>ob s (t) (solid line), with the best- 
fitting function for the least-damped first-order QNM, 4>Bt (t) 
[Eq. (|5ip ] (dotted line). The least-damped QN frequency is 
so = —0.50 — 2.18i [see Eq. (|49}]. The observer's position is 
Zobs = 5.0. 



be composed of three types of components: the second- 
order QNMs, second-order power-law tail, and first-order 

(2) 



order QNMs. To evaluate the excitation coefficient Co 
of the least-damped mode, we adopt a method used in 
Ref. [25[: first, we calculate a quantity 



we introduce a fitting function 
Eqs. ||3SJ| and (03J)] 



(t) defined by [cf. 



£(1) 



^obs 



(*) 



dt. 



(50) 



where t\ and t 2 are fixed lower and upper limits of in- 
tegration, and cf>g) (t) is a fitting function defined by [cf. 
Eq. EZJ] 



b {2) 
"fit 



4V(t) = Coe 80 ^"^) + C* e s S(*-*ote) 

= 2|C |e^(*-^) cos[O (i - x ohs ) - arg(Co)], 

(51) 

where Co = ICole 1 " 8 ^ ' is the excitation coefficient of 
the least-damped first-order QNM, and sq is the fre- 
quency of the least-damped first-order QNM for the 
Poschl- Teller potential [Eq. (|49|) ]. The fitting parame- 
ter is Co (complex). Then, we search for the value of 
Co = Co,fit that minimizes . The value of t\ is taken 
to be larger than x bs, because the overtone QNMs will 
not be negligible at very early times. We chose t\ = 6.0 
and t-2 = 25.0, obtaining C 0)fit = -0.176 - 0.22H. For 

comparison, we plotted in Fig.[3]the function (f$ (t) with 
the best- fitting parameter Co = Co, fit- It is clear that the 
least-damped mode dominates the numerical waveform 
except for very early times 4.0 <, t <> 5.5, for which the 
overtone modes will not be negligible. 

In Fig. 0Ja), we show the observed waveform of the 

second-order perturbation, <f>^ s (t) = $ 2 \t, x Q h s )- From 

the result of Sec. Ill, we expect that (p^ a (t) should 



\ 2Q (t) = a 2Q C m Me 2s ^- x °^ 

+ |a 2 Q|^oo,fite 27o(t -^ ) 
i2) (t) = a T H &t ( 



c.c. 



1 



1 



y fit,XV 



v t ^obs ^2 t ~\~ 3^obs ~\~ ^ 



(52) 

(53) 
(54) 

'(55) 



where Coo, fit, £>oo,fit, and -f/fit are the excitation coeffi- 
cients of the second-order QNMs and the power-law tail 
defined by Eqs. (gDJ), (gl|), and @J with C = C ,fi t 
and C„ = for n > 1. The fitting parameters are 

( 2) 

Cq (complex) and two scaling parameters a 2 Q (com- 
plex) and (Xt (real). The scaling parameters a 2 Q and 
cut are introduced to evaluate how accurately the asymp- 
totic approximation predicts the amplitude of the second- 
order components. Searching for a set of the best-fitting 



parameters 



( 2) 

,o.t,Cq )fit as done for 



6 (1) 

y obs> 



we ob- 



tain a 2Qfit = -0.0277 + 0.21H, a T M = 0.528, and 

C (2) Rt = -0.00799 - 0.0102i. We find an excellent agree- 

ment between the best-fitting function <p^ t (t) and the 

numerical waveform (f>^ s (t), as shown in Fig. SJa). The 
fact that |a2Q,fit|, o;T,fit < 1 implies that the asymptotic 
approximation overestimates the true amplitude of the 
evolved perturbation by a factor. This is because the 
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asymptotic approximation breaks down near the poten- 
tial peak, \x\, \x'\ < 1, where the source term S^ 2 ' gives a 
large contribution to <f>^ [recall our choice of F(x) and 
see Eq. ([501) ]. It should be noted that the asymptotic ap- 
proximation correctly predicts what types of components 
construct the evolved second-order perturbation, which 
we really want to know in this study. 

As explained in Sec. Ill, what we are interested in is the 
second-order waveform from which the first-order QNM 

(2) (2) 

component is removed, o b B (*) — <%t iq(*)' We plot this in 
Fig. HJb) , comparing with the "truly second-order" COm- 
^l (21 

ponents, (p^ t 2 q W + ^ t T (f). We find that the tail com- 
ponent dominates the extracted waveform at all times, 
but nevertheless the second-order QNM component is not 
negligible at t ^ 9. In Fig. BJc), we plot the numerical 
waveform with all the components other than the second- 

(2) (2) C2) 

order components removed, </> obs (i) — [(f> Rt T {t)+(jr Rt ^(t)}. 
We find that the extracted numerical waveform can be 

(2) 

well reproduced by a superposition 0g t 2 q(0 01 tne least- 
damped second-order QNMs (i.e., sum tone oc e 2s °* and 
difference tone oc e 2lat ). 



C. Result 2: The fully nonlinear field 

In Fig. [5l we compare the result of the fully nonlin- 
ear calculation of Eq. ([7]), $ o bs(0 = £obs), with the 
numerical waveform of the first- and second-order pertur- 
bations, </>ob s (i) an d ^obsW- At early times, t <, 12, no 
significant difference is seen between 3> bs(i) and 0obs W- 
This means that the first-order theory is sufficiently ac- 
curate and thus any higher-order correction is unimpor- 
tant in this stage. At late times t 12, however, <^ ob 
fails to predict the behavior of ^obs- Surprisingly, it is 

(2) 

the second-order perturbation obs that shows a good 
agreement with ^obs in this stage. This means that the 
second-order perturbation dominates the late-time be- 
havior of the nonlinear waveform, and the perturbations 
of all the other orders are not significant. We have also 

(2) 

confirmed that </> obs continues to agree with $ Q b s as late 
as t = 50, which strongly suggests that third and higher- 
order perturbations never surpass the second-order one 
in amplitude for all late times. 



V. CONCLUSION AND DISCUSSIONS 

In this study, we have investigated the nonlinear evo- 
lution of black hole ringdown using second- and higher- 
order perturbation theory and a simplified nonlinear field 
model for black hole metric perturbations. We have 
proven that second-order QNMs. whose existence has 
been predicted by recent works [15|, [l6[, do appear in 
the evolved second-order perturbation. As a bonus, we 
have discovered that the second-order QNMs are accom- 
panied by a new type of power-law tail. This power-law 
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FIG. 4: Comparison of the numerical waveform and the 
best-fitting functions for the second-order perturbation (jr- 2 ' 
at x hs = 5.0. (a) The numerical waveform 4> b s (t) (solid) 
and the best-fitting function for the least-damped first-order 
QNM <f>$ <1Q {t) [Eq. (dotted line), (b) The numeri- 

cal waveform with the first-order QNM removed, 4>obs(t) ~ 
^'^(f) (solid line), compared with the sum of the best- 
fitting functions for the second-order QNMs and power-law 
tail ^ (2) ^ -l ^ (2) 



"fit,2Q 



(t) + <j>£ t > T (t) [Eqs. flS3} and JMJl] (dotted line), 



and the best-fitting function for the second-order power-law 



tail only, 



i(2) 



(t) (dashed line), (c) The numerical wave- 



form with all the components other than the second-order 

QNMs removed, ^(t)-[<t>^ tT (t)+^ 1Q (t)] ( solid line )> and 
the best-fitting function for the second-order least-damped 
QNMs, (j>fit2Q(t) (dotted line), which is a superposition of the 
"sum tone" mode oc e 2s<)t (dashed line) and "difference tone" 
mode oc e 2jot (dot-dashed line), where so = 70 — iQo is the 
frequency of the least-damped first-order QNM [see Eq. (|49p ] . 
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t 



FIG. 5: Numerical waveforms of the fully nonlinear field 
(solid line), the first-order perturbation (dotted line), and 
the second-order perturbation (dashed line). Also plotted 
are the first-order perturbation with the second-order correc- 
tion added (dot-dashed line) and the best-fitting function for 
the second-order least-damped QNMs obtained in Sec. IV. B 
[Eq. (J53J)] (dot-dot-dashed line). 

tail, to which we have referred as the second-order power- 
law tail, decays more slowly than any of the tails in the 
first-order theory, and even dominates the fully nonlin- 
ear evolution at late times. In other words, the first-order 
(i.e., linear) perturbation theory fails to predict the late- 
time evolution of the ringdown, for which higher-order 
corrections must be taken into account. Also, we have 
shown that the behavior of the second-order tail is de- 
termined by the asymptotic form of an effective source 
term, i.e., a term expressing the nonlinearity of the sys- 
tem. Since the asymptotic form of the source term in our 
model is the same as that in black hole perturbations, the 
second-order power-law tail as well as the second-order 
QNMs will certainly appear in real black hole ringdown. 
These nonlinear components could open a new precision 
science in gravitational-wave studies. 

Interestingly, the surprising failure of linear pertur- 
bation theory has been reported just recently for other 
nonlinear fields [2(| [U Ha H3| • Power-law tails of non- 
linear orig in have been first discovered in the Skyrme 
model [26], and later in a spherically symmetric Yang- 
Mills field on flat and Schwarzschild backgrounds [27| (in 
fact, their nonlinear tails are of third-order perturbations, 
because the second-order ones cancel out; as mentioned 
in Ref. these seem to be special cases). In both 

cases, the nonlinear tails decay more slowly than the lin- 
ear ones, and do dominate the fully nonlinear evolution at 
late times. From this fact, we suspect that our second- 
order power-law tail would have close relation to these 
nonlinear tails. However, we emphasize that the nonlin- 
ear tail of metric perturbations has been for the first time 
discovered in this study. 

In our future works, we plan to address the following 
open questions: 

(i) Do the second-order QNMs and power-law tail ap- 



pear in observable quantities, i.e., in gravitational 
waves? As explained in Sec. II, the perturbations 
<fP^ and </>( 2 ) of our model field correspond, re- 
spectively, to the master functions ip^m an d 
of the perturbations of the Schwarzschild geom- 
etry [see Eqs. ©, ©, ©, and (fl0j)]. There- 
fore, it is straightforward to apply our analysis 
used in this study to the evolution of the mas- 
ter functions, although the computations will be 
much more complicated. It should be noted, how- 
ever, that these master functions themselves arc 
not observable quantities. In order to convert them 
into truly observable quantities, we have to perform 
a gauge transformation from the Regge- Wheeler 
(RW) gauge to the asymptotically flat (AF) gauge, 
and then compute the transverse-traceless (TT) 
parts of the metric perturbations in the AF gauge 
JT5L [Hj]. We will have to prove that the second- 
order QNMs and power-law tail remain to appear 
in the TT parts of the metric perturbation. Giving 
a full-nonlinear proof should be a nice challenge for 
numerical relativity. 

(ii) Do our results remain unchanged under more re- 
alistic initial conditions? In the present study, we 
have assumed initial perturbations localized around 
the potential peak. In order to deal with astrophys- 
ically important problems, we need to extend our 
analysis to more general initial conditions. One of 
the most tractable problems among them will be 
the head-on collision of equal-mass Schwarzschild 
black holes [3(J, [M| . The evolution of the second- 
order perturbations during the collision was already 
studied [2(| using Misner's initial data [32|], but the 
second-order QNMs and power-law tail were not re- 
ported in that study. Therefore, we must discuss 
this issue using the same initial data. We expect 
that the effect of initial perturbations located far 
away from the potential peak would not be signif- 
icant, since the effective source term is suppressed 
at r ^ 2M by a factor of oc r~ 2 . 

(iii) Is it possible to detect the second- order power-law 
tail ( if it exists ) with future gravitational-wave de- 
tectors? To address this issue, we need to estimate 
the energy carried out by the second-order tail as- 
suming astrophysically realistic events (e.g., binary 
black hole mergers), as already done for second- 
order QNMs [la, LL6|]. Also, we have to develop a 
data analysis to extract the second-order tail from 
detector outputs. While the matched filtering tech- 
nique has been developed for extracting QNMs [H[ , 
there seems to be only a few techniques for power- 
law tails [HI]. If the second-order tail is detectable, 
it could be used to distinguish true signals and spu- 
rious ones in black hole ringdown search, where 
fake reduction and event identification are crucial 
Bill- 
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APPENDIX A: CONSTRUCTION OF GREEN'S 
FUNCTION 

In this appendix, we summarize the derivation of the 
flat part Gf and QNM part Gq of the Green's function 
using the asymptotic approximation [T^ |. 

First, we take the Laplace transform of G(r; x, x') with 
respect to r: 



G(s; x, x') 



r G(T;x,x') 



(Al) 



This is an analytic function of s defined for Re(s) > 0. 
Taking the Laplace transform of Eq. (fT8"|) as well, we 
obtain 



l-s 2 



V{x) G(s;x,x') =8(x-x'). 



(A2) 



To construct the solution G, we introduce two homoge- 
neous solutions to Eq. (|A2|) . yi,(s,x) and y R (s,r), which 
satisfy the boundary conditions 



Vl{s,x) 
Vr{s, x) 



A L {s)e sx + B L (s)e- 



A R (s)e~ 



B R {s)e s 



(x — > — oo) 
(x — > +oo) 

(x — ► — oo) 
(x — ► +oo) 



(A3) 
(A4) 



where Al, Bl, A r , and Br are functions of s deter- 
mined by the form of the potential barrier. Physically, 
\B L (R)(s)/A L ( R) (s)\ and \l/A L(R) (s)\ mean, respectively, 
the reflection amplitude and the transmission amplitude 
of a wave incident from x — > +oo(— oo) with unit ampli- 
tude and a frequency lo — is. With these homogeneous 
solutions, G(s; x, x') is constructed to be 



1 



G(s; x, x') 



-yL(s,x)y R (s,x'), (x<x') 



W [°> (A5) 
^77^yyi(s, x')y R (s, x), (x' < x) 



where W(s) is the Wronskian defined by 

W(s) = y L (s, x)d x y R (s, x) - y R (s, x)d x y L (s, x). (A6) 

It can be easily shown that W(s) is independent of x. 
Evaluating the right-hand side of Eq. (|A6|) at x — > ±oo 
gives W(s) = -2sA L {s) = -2sA R (s). Equation (|X5]l 




-r- Res 



FIG. 6: Schematic illustration of the contour for the inte- 
gration in Eq. (|A7|I . The solid line represents the original 
contour. The zigzag line is a (possible) branch cut for the 
integrand. 



is the only possible choice of G for the boundedness of 
time-domain solutions p], [ll| . 

Taking the inverse Laplace transform of G, we have 



r e + io ° ds 

G(t;x,x') = / e sr G(s; x, x') 

' ioo 2tti 

'■ e+ ' io ° dse ST 



Je-ioo 4nisA L (s) 

x iyL(s,x)y R (s,x'), (x<x') 
\yL{s,x')y R (s,x), (x'<x) 

where e > 0. This integral can be evaluated by closing 
the contour of integration with an infinite semicircle in 
either the left or right half in the complex s-planc (if G 
has a branch point at s = 0, which is true in many im- 
portant cases, the infinite path is taken to circumvent 
the branch cut; see Fig. 0. For this purpose, we ana- 
lytically continue the integrand in Eq. (|A7|) into the left 
half of the complex s-plane, Re(s) < 0. The analytic 
continuation of the integrand has an infinite number of 
poles {s n }, which can be shown to be equal to the ze- 
ros of Al(s). This implies that analytic continuation of 
the homogeneous solution yL(s,x) satisfies the "outgo- 
ing boundary conditions" (|26[) at s = s n , meaning that 
y n (x) = yL,(sn,x) are the first-order QNMs. 

It is generally hopeless to carry out the integration in 
Eq. (|A7[) exactly since neither y^, y R , nor Al can be 
written as a simple function of s. Nevertheless, it is pos- 
sible to carry out the integration approximately by using 
the fact that yL and y R approach to exponential forms 
far away from the potential, |x| 3> 1 [see Eqs. (|A3[) and 
(|A4[) ]. We call this "asymptotic approximation" after 
Andersson 0. From Eqs. (fATij) and (fX4| . Eq. (|a7|) is 
approximately written as 
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G(r; x, x) 



e+zoo 



ds 



s(t-(s-x')) 



(a:' < -1) 



+ A L (sf ' 1 ^ J 



(A8) 



where we have assumed that x ^> 1 for all the cases. 

Now we choose the contour of integration according to the convergence of the integrand in Eq. (|A8j> at \s\ — > oo. 
For x' <C — 1, the integrand converges in the left half-plane if r — (a; — x') > 0, and in the right half-plane otherwise. 
Taking a contour in the half-plane where the integrand converges and applying the residue theorem, we obtain 



G(t; x, x') 



E 



2 S,vAr (S n ) 

n— ' v v y 



o. 



(r > a; — a;') 
(t < x — x') 



(A9) 



where A' L (s n ) = (IAl(s) / 'ds\ s — Sn . For x 3> 1, the integrand involves two exponentials with different arguments, and 
thus a preferred choice of the contour is not obvious. Following Ref. [13], we choose to close the contour in the left 
half-plane if r — (x + x') > 0, and in the left half-plane otherwise. This yields 



G(t;x,x') « < 



2 s„j4'r(s n ) 

n— ^ 

1 

2' 



0. 



(t < x + x' and r > |a; — x'\ 
(t <\x- x'\) 



(A10) 



Combining Eqs. (|A"9j) and (|A10|) together with Eq. ((2B]| gives G F and Gq in Eqs. (JM]) and (|25]l. It is noted that the 
"tail part" Gb can not derived from the asymptotic approximation, since the asymptotic approximation breaks down 
in the low- frequency region (\s\ rs 0). 
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